Dynamics of spherically symmetric spacetimes: Hydrodynamics and Radiation 
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Abstract 



qq , Using the 3+1 formalism of general relativity we obtain the equations governing the dynamics 

of spherically symmetric spacetimes with arbitrary sources. We then specialize for the case of 
perfect fluids accompanied by a flow of interacting massless or massive particles (e.g. neutrinos) 

\£> . 

, which are described in terms of relativistic transport theory. We focus in three types of 

coordinates: 1) isotropic gauge and maximal slicing, 2) radial gauge and polar slicing, and 3) 



isotropic gauge and polar slicing. 
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I. INTRODUCTION 



One of the most fascinating phenomena in gravitational physics is that of gravitational col- 
lapse. Notably, the gravitational collapse of astrophysical bodies culminating in the formation 
of black holes. 

The historical controversy on the final fate of gravitational collapse of compact objects 
such as white dwarfs and neutron stars raised by the discovery of maximum mass limits and 
the subsequent stability analysis led to find, at least for the most ideal configurations, definite 
answers and concrete predictions depending on the initial conditions and the equation of state 
of matter (see Refs. for a review). The simplest situation describing the gravitational 
collapse ending in a black hole formation is that of a spherical ball of pressureless and ho- 
mogeneous fluid (the well known Oppenheimer-Snyder dust collapse ||]). The solution shows 
the appearance of an event horizon revealing thus the formation of a black hole after a finite 
proper time. 

Since that pioneering investigation, much has been done for more complicate and realistic 
initial configurations. On one hand, the accelerated development in the area of computation 
and the recent advances on the numerical analysis of Einstein equations have made possible 
computing the dynamics of rather complex spacetimes faster and for longer term evolutions. 
On the other hand, the advances in particle and nuclear physics have led to a better knowledge 
of the conditions of matter at high densities, providing then more realistic models for the matter 
in cores and neutron stars. 

Most of the recent analysis of gravitational collapse leading to a black hole formation have 
been done in spherical symmetry and following a more "modern" point view in light of the 3+1 
formulation of general relativity and other formulations better adapted to numerical stability 
(see Ref. and the references therein). Among these investigations, there is the one of 
Shapiro & Teukolsky || who studied the collapse of polytropes by imposing an isotropic gauge 
and maximal slicing coordinates, with initial conditions provided by the Tolman-Oppenheimer- 
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Volkoff equation of hydrostatic equilibrium. Later, a similar study, was performed by Schinder, 
Bludman and Piran 0] in comoving coordinates with a polar slicing, and by Gourgoulhon |8|,^] 
in radial gauge and polar slicing, this latter improving previous analysis by the incorporation 
of realistic equations of state for the nuclear matter. 

The effects of the coordinate choice and the slicing condition on the time evolution has 
been exhibited for two of the most popular choices (isotropic and radial gauges with maximal 
and polar slicings) by comparing with the analytical solution of Oppenheimer-Snyder (OS) 
For the time being, the discussion of gauges is postponed to Sec. VII. 



The analysis of gravitational collapse of compact stars and iron cores would not be complete 
if the influence of neutrinos were not taken into account. This has been recognized by a long 



list of authors since the precursor investigations of Colgate and White [12], and May and 



White [O], who analyzed the effect of neutrinos in supernovae explosions. Later Wilson 14], 
performed full general relativity computations with a neutrino flow described in terms of 
the relativistic Boltzmann equation. Wilson's analysis included electron and muon massless 
neutrinos assuming that the corresponding antineutrinos contributed in the same basis. The 
interaction of neutrinos with the star's fluid was described by an opacity function. Unlike 
previous studies, Wilson's found that the heat conduction by neutrinos is not sufficient to 
eject any material from a collapsing star. In all those studies black hole formation were not 
analyzed but only evolution configurations terminating in stable states corresponding to white 
dwarfs or neutron stars. 



An updated analysis were carried out by Mayle, Wilson and Schramm [15] using a Boltz- 



mann code and for a large set of mass configurations. Neutrino signals from various species 
were analyzed within time scales of ~ 1 s after the supernova explosion. It is to be men- 
tioned that previously Saenz and Shapiro fl(| had computed a non-spherical Quasi-Newtonian 
collapse acompanied by neutrino and gravitational radiation. 

Burrows and coworkers have also analyzed during the last twenty years the mechanism 
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of Type II supernovae (SNII) explosions and the role of neutrinos (see [ p^?|JT8| and references 
therein). Among these investigations, we find an interesting model of long term neutrino 
emission from the hot protoneutron star phase to the final outcome of a stable cold neutron 



star 19 . 



Many other recent investigations have confirmed and improved in several aspects previous 



findings on SNII (see [20-25] and references therein). For the case of a core collapse leading 
to a black hole two scenarios are recognized p2| . One called early black hole formation which 
is generically associated with accreting protoneutron stars which form from the collapse of 



degenerate cores of massive stars 25,2™. The accretion of some tenths to one solar mass can 
last a second, and the exceeding of the maximum mass drives the protoneutron star into a 
black hole collapse in a typical time scale of ~ 0.5 ms. In this scenario the neutrino signal 
is abruptly cut off after the black hole forms, and the typical neutrino luminosities prior the 
cutoff are ~ 10 52 erg/s per flavor. 

The second scenario called late black hole formation typically arises by a softening of the 
high-density equation of state of the protoneutron star |27-^0|. The phase transition from the 



neutron star matter to a more exotic state which include kaon condensates [31,32], or hyperon 



condensates [33| can lower the maximum allowable mass to ~ 1.5M 28. 3C], driving thus a 
stable protoneutron star to an unstable regime and finally to a collapse into a black hole. This 
kind of core collapse can last ~ 10 s before the cutoff and the luminosity of neutrinos is ten 
times lower than the luminosity of the early case. 

It is encouraging to note that for a SNII at a distance of ~ 10 kpc which explodes within 
the early scenario, SuperKamiokande can probe v e masses down to 1.8 eV by comparing the 
arrival times of high and low energy neutrinos within the reaction v e +p — * e + +n in a Cerenkov 
detector (see Ref. [^] for details). 

In fact the very recent announcement on the measurements of solar neutrinos from the 
decay of 8 B by the Sudbury Neutrino Observatory (SNO) |34j via charged current interactions 
and by the elastic scattering of electrons reveals that neutrinos could be changing flavor as 



they travel from the sources to the Earth. This discovery if confirmed could corroborate the 
oscillating behavior of neutrinos and therefore their massive nature. The fluxes measured of 
the different flavors are in close agreement with the predictions of the solar models. The SNO 
experiment then implies that the upper limit of the mass squared difference between the v e 



and the v T or is less than 10 eV [34|. This result when combined with the current bounds 
on m Ue of 2.8 eV and Am^ „ t (assuming neutrino oscillations) provides a limit for the sum of 
the masses of the three neutrino species in the range [0.05,8.4] eV plf . 

One proposal to measure the v T and Vn masses indirectly and that can corroborate the 



SNO findings is the one which uses a time-of- flight technique [22], for neutrinos emitted in the 
early black hole formation scenario discussed above. The point is that if neutrinos are massive 
then there is a delay (relative to a massless neutrino) in the cutoff of the neutrino signal as 
measured on Earth after the black hole forms, and it is given by At ~ {m u /E v ) 2 for distances 
of ~ 10 kpc. This delay can affect the event rate measured in a detector. The conclusion is 
that assuming luminosities L ~ 10 52 erg/s per flavor at the cutoff time, SuperKamiokande can 
probe e-neutrino masses as small as 1.8 eV for T Ue ~ 3.5 MeV, whereas the OMNIS |35|] or 
SNO detector can detect m VlMjVT masses as small as 6 eV for T v „ T ~ 8 MeV. 

A collapse scenario which is rather different from the above consists in the accretion of 
matter by an old neutron star near the maximum mass limit. Gourgoulhon and Hansel [Q 
have analyzed the neutrino emission during the collapse to a black hole within this scenario via 
nonequilibrium j3 processes, assuming that the nuclear matter is transparent for neutrinos (i.e., 
opacities were neglected). Instead of using neutrino transport, a regularized geometrical-optics 
model adapted to massless neutrinos was adopted. This model was thus intended to provide 
upper bounds in the neutrino burst. The collapse lasts typically a millisecond (the time that 
takes place the black hole formation), and in the most favorable conditions the total energy 
of v e and antineutrinos is ~ 10 51 erg, while the energy of the corresponding neutrinos is 
several orders of magnitude lower. This is even lower for the v T and v T neutrinos. The main 
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conclusion is that a collapse of this kind at a distance of ~ 10 kpc would be undetectable by 
the current neutrino detectors. 



Finally, another scenario which has been analyzed in the past is the dynamics of collisionless 
gas of particles which mimic spherical star clusters, and the possibility of a cluster collapse into 



a supermassive black holes [38-jilq|. The motivation was to provide a theoretical description for 
the formation of supermassive black holes that could exist in the centers of galaxies. Recently, 
a similar study which include gamma-ray bursts, is the one analyzed by Linke et al. |43[| , where 
the collapse of supermassive stars (M ~ 1O 5 M - 1O 9 M ) with emission of thermal neutrinos 
is considered. In that work, the spacetime is foliated by outgoing null hypersurfaces rather 
than using a 3+1 foliation of spacetime. 

In view of the different scenarios of gravitational collapse available today and the miscella- 
neous predictions within each of them it is worth pursuing the investigations along these lines. 
Only in this way there will be at hand a large set of models which the forthcoming (e.g. |35| ) 



and recent observations [34] will validate or rule out. 

Although the paper is written in the same spirit of various papers which deal with the sys- 
tem of equations Einstein-Hydrodynamics-Boltzmann, several aspects distinguish from them. 
For instance, most formalisms treat neutrinos as massless particles, except perhaps the one of 



Harleston and collaborators [ 44 , 45 1 . Here massive particles are considered from the onset and 
previous equations are recovered in the massless limit. The relativistic Boltzmann equation is 
written in terms of 3+1 variables for generic spacetimes. This had done in the past only in 
spherical symmetry. Therefore, the relativistic Boltzmann equation presented here is coupled 
from the onset to the 3+1 Einstein's equations. That is, the curvature effects appear in terms 
of the lapse, the shift, the extrinsic curvature and the three-metric. The hydrodynamic equa- 
tions are derived also in the context of the 3+1 formalism and they couple to the neutrinos 
via collision integrals. In particular the equation for the velocity field of the fluid is written 
in several forms each of one is useful whether one uses different numerical methods. At this 
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regard, a general relativistic Euler equation is presented using the tetrad formalism. Its quasi- 
Newtonian form allows an easy interpretation of several terms, and reduces to well known 
equations in various limits. Such an Euler equation turns to be better adapted for spectral 
methods than the equation for the momentum current density |8],|[J . The system of equations 
are then specialized for spherical symmetry and written in three different coordinates: 1) 
isotropic gauge and maximal slicing, 2) radial gauge and polar slicing, and 3) isotropic gauge 
and polar slicing. 

This is the first of a series of papers where the gravitational collapse of various kinds of 
matter will be analyzed. 

The paper is organized as follows: in Sec. II we present succinctly the 3+1 formalism 
of general relativity rather more to fix the notation than to give a detailed description. In 
Sec. Ill we consider the case of two interacting sources of matter: a perfect fluid and a flow 
of relativistic particles described in terms of relativistic transport theory. Section IV treats 
the relativistic transport theory. Section V deals with therodynamics. In section VI spherical 
symmetry is considered and in section VII three coordinate choices and slicing conditions are 
analyzed and discussed in light of the previous studies. Finally we conclude with some remarks 
and the plans for the forthcoming investigations along this line. 



II. THE 3+1 FORMALISM OF GENERAL RELATIVITY 

One of the most popular reformulations of general relativity when tackling numerical prob- 
lems is the 3+1 or Adison-Deser-Misner (ADM) formulation. We shall not enter into the details 
of the derivation of the equations (see Refs. p6|-[50|]) but rather discuss the general idea in order 
to fix the notations. 

The main idea is as follows: under general assumptions (see |47 ,|5C|l and reference therein 
for details) a globally hyperbolic spacetime (M 4 ,^^,) can be foliated by a family of space- 
like hypersurfaces S t (Cauchy surfaces). Each hypersurface represents a Riemannian sub- 
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manifold (M 3 , h^ ) endowed by an induced metric (the 3-metric) . It is then assumed a local 
coordinate system (t, x l ) for the spacetime, the spatial part (x % ) represents a local coordinate 
system for Ej, while t is a global time function that parametrizes E^. The embedding of St in 
spacetime is completed by the extrinsic curvature of E^. This is defined by 

:= ——Cnhpu , 

where Cn stands for the Lie derivative along the normal to Ej and fy^, := g^ u + n^n v . 
The vector field is time-like (n^n^ = — 1) and the convention used for its components with 
respect to the coordinate base adapted to the spacetime foliation is as follows: 

This convention means that points towards the future. Since is a unit time-like vector, 
it is customary to interpret as the four- velocity of the so-called Eulerian observer Oe- The 
scalar quantity N (the lapse function) represents thus the rate at which Oe sees the flow of 
its proper-time as compared with the intervals between two neighboring hypersurfaces Ej and 
Et+dt- The 3- vector N l (the shift-vector), represents the coordinate 3- velocity at which the 
Eulerian observer moves with respect to the coordinates (t,x l ). In this way, the four-metric 
reads 

ds 2 = -(N 2 - WN^dt 2 - 2N l dtdx i + h ij dx i dx j . 

We emphasize that in many studies the convention on the sign of the shift vector is different 
from the one adopted here. This has to be taken into account particularly for purposes of 
comparison between the final form of equations in spherical symmetry presented here and the 
corresponding equations of the references which employs the opposite sign. 
A useful formula for Kij obtained from Eq. (|l]) is 

Kij = -V inj = -m% = -±(^> + 3 VjNi + s ViN^ . 

where 3 Vj stands for the covariant derivative compatible with hij. This is to be regarded as 
an evolution equation for the three- metric hij. 



The trace of the extrinsic curvature is simply given by, 



K := -V a n a 



Another useful quantity is the acceleration of Oe given by 

■= n u V v n» = 3 V^[lniV] , 



which allows for the lapse the interpretation of the acceleration potential for 0e [ 49 1 • 

The orthogonal decomposition of the energy-momentum tensor in components tangent and 
orthogonal to T, t leads to 



The tensor S^ v is symmetric and often called the tensor of constraints] is the momentum 
density vector and E is the total energy density measured by the Eulerian observer 0e- Both 
and J M are orthogonal to n/ 4 . 

We emphasize that for the specific applications we will study, T^ v will be the total energy- 
momentum tensor of matter which can be composed by the contribution of different types of 
sources: 

i 

This means that 

i i i 

The projection of Einstein equations = 4irGo (2T^ V — T a a g^ u ) in the directions tangent 
and orthogonal to E$, followed by the use of the Gauss-Codazzi-Mainardi equations leads to 
the 3+1 form of Einstein equations: 

3 R + K 2 - KijK lj = 16ttG E , 

known as the Hamiltonian constraint. 



3 V t K l t - 3 V t K = SvrGo J { , 

known as the momentum constraint equations. 
Finally, the dynamic Einstein equations read 

d t K i j + N'-diK^ + K\djN l - K l foN* + 3 V* 3 VjN - 3 R i j N - NKK i j 



AnGnN 



(S - E)5 i j - 2S) 



where S = S\ is the trace of the tensor of constraints, and all the quantities written with a '3' 
index refer to those computed with the three-metric h™. Moreover, under the 3+1 formalism 
tensor quantities tangent to T,t use the three-metric to raise and lower their spatial indices. 
The equations @ and (|l^) are the set of the Cauchy-initial-data evolution equations for the 
gravitational field subject to the constraints Eqs. ([!(]) and 

An evolution equation for the trace K is obtained by taking the trace in Eq. (12): 



d t K + N l d t K + 3 AiV - iV ( 3 R + if 2 ) = AttG N [S - 3E] . 

where 3 A stands for the Laplacian operator compatible with h{j. 
This can be simplified by using Eq.([io|) to give 



d t K + N l diK + 3 AN - NKijK ij = 4ttG N [S + E] . 
To complete the system of equations, we consider also the matter equations 

V M T^ = , 

which according with Eq. (||), these can be written as 

where T^ u is the energy-momentum tensor of certain fields that are collectively labeled by ip 
and T v := V^T^ are the "forces" exerted by the external fields (fields other than ip). For 
instance, we shall consider the case where the total energy-momentum tensor is given by a 
combination of a perfect-fluid and a radiated flow of particles: 
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so that T^ V = Tpp and T^ t = thus, the energy-momentum tensor of the perfect-fluid 
alone will not conserve by separate; J- v will represent the "forces" of the radiated flow acting 
on the perfect fluid in form of collisions. 

The energy- momentum conservation equations ( |l6| ) can be written in 3+1 form as well. 
The projection of Eq. (|l6|) along leads to the energy conservation equation, 



8 t E^ + N 1 ^ + — 3 Vi {N 2 J l A = N (S l jK i:j + E^K) + Nn v T v 



Explicitly n v T v = -NT 1 . 

On the other hand, the projection of Eq. ( [1~6| ) on leads to the momentum conservation 
equation, 

dtjf + N l d t jf + jfdiN 1 + N( 3 V^S\) = NKjf - + E^S\) 3 V t N - . 

where 



'T = hipf* = -NiF + hyp . 



A. Tetrads 

In many applications the use of tetrad components of tensors (hereafter physical compo- 
nents) are better adapted to a problem than the coordinate components. This will be the case 
when writing the equations of relativistic transport for the radiated flow and the equations of 
motion for the matter. 

In the context of the 3+1 formalism the tetrad we use is the local tetrad of the Eulerian 
observer which is given by {n^e^}, that is, by the time-like vector normal to Et and by a 
triad on E^. In covariant notation that tetrad is given by: 

_ u 9 
11 



where d/dx v denotes the coordinate basis of the spacetime, and q^s are the tetrad coefficients 
that allow the normalization. For instance, it turns that ej^ = n^. 



The inverse relationship of Eq. (|2l]) is given by 



dx^ 

where the coefficients effl are related to by the completeness relations ei"^^ = ^fy, 
and ei Q) ^ a) = 6% . 

A tetrad is not uniquely defined. The Lorentz invar iance SO (3,1) leaves the freedom on 
the choice of the six parameters that rotate and boost frames. As mentioned, a convenient 
choice is as follows, 



e? = N 



'hj — e « e j 



The inverse relations are, 

i _ N l 

io = , 

it ) = N- 1 , 

here ef^ = 5^ e\ l) q l {j) = 5 (l) {j) . 

This choice of a tetrad is compatible with the 3+1 decomposition of the four- metric (||), 
so that 

with given by Eqs. (p3|) — (26) and 7/(a)(/9) stands for the Minkowski metric. 
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Trivial Lorentz transformations of the chosen tetrad relates the Eulerian observer with 
other possible frames. Moreover, the best choice for a triad ep on £j will depend on the 
particular coordinates used on the hypersurface E< (see Sec. VI for the case of spherical 
symmetry) . 

Finally, the transformation law for the components of tensors tangent to St from the 
coordinate base to the triad is as follows, 





= efV , 


J(i) 


= l{i) J i i 


(?) 


= e^S* 


q(i) 


= e? ^ 



The inverse relationships are obtained from above in the obvious way. The triad indices [i.e., 
spatial indices within '()'] are raised and lowered with <5 9-s (i.e., the triad-covariant and triad- 
contravariant components of 3-tensors are identical to each other). Four-tensor components 
transform in a similar way using the four-dimensional tetrad coefficients and V(p)(i/) (resp. 
^MM) to lower (resp. raise) indices. 

The use of the tetrad formalism will be useful to recast the 3+1 matter equations and later 
on to write the relativistic Boltzmann equation in a useful manner. For instance Eq. (|l£ 
reads, 



d (t) E^ + 3 V (l - } Jf - E^K + 2JW 3 V {i)1 y - j,S®WK, 



W(i) 



where du\ = n^d^. 

On the other hand, the momentum conservation equation (|19|) can be written as 



W(j) 



+ 



V (i)^ 



o 



(i) 

m) 



K 



(0 
(0 



. 3 



We remind that the covariant-derivative components with respect to a tetrad uses the 
four-Ricci rotation coefficients (RRC) as a connection. We define them as follows, 



13 




The three-RRC have an identical expression by restricting the above definition to pure spa- 
tial indices and using the three-covariant derivative. While the above definition requires the 
Christoffel symbols, these can be avoided by using the representation of the RRC in terms of 
the structure constants 

Although the 3+1 Einstein equations can be written following a tetrad approach [^], for 
our purposes this will not be necessary and thus we will not pursue the issue here. 
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III. PERFECT FLUIDS WITH SOURCES 

For the specific applications we have in mind, a combination of a perfect-fluid and a 
radiated flow will be considered. Then in Eq. ( |l7|) we assume 

T^ = (p + p)u^+pg^ . 

For the moment the form of the radiating part T^ u is not specified. This will be treated in 
detail in Sec. IV. 

The corresponding 3+1 matter variables of the fluid are, 



Epf 


= (p + p)T-p , 


r(0 
J PF 


= (E PF +p) 3 U^ , 


aim 

J PF 


= {E PF + p) 3 U il) 3 U^ + $WCj) 


SpF 


= (£ PF +p)( 3 [/«) 2 + 3p , 


r 


, = _ n ^ = u {t) = h_ (3 [/ (i))2 



where 

3 C/ « . = = _L_ (V _ N i) = 1 ( V (i) _ N (D) 

r n v j n v j 

V 1 := u l /u l . 
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The equation for conservation of energy (p6|) applied to a perfect fluid with sources then 



reads @|||-|| 



d (t) E + 3 v (0 



(E + p) 3 U®} +{E + p) [2 3 U {j) a {j) - 3 C/« 3 U^K t 



K 



On the other hand, the momentum conservation equation (37) applied to a perfect-fluid 
with sources can be written as an Euler equation for the velocity field. This reads |48|,52|-j54|], 

1 



d (t) 3 c/ (i) + 3 u ij) 3 v m 3 u^ 



U) 



Ep F +p 



l d^p+ 3 U^d (t) p 



Ki) 



+ 3 U^ 3 U«Ua (l) - 3 U^K m ^ ' 



where 



(*)(0 



K 



+ 1 f SjjW^t) _ 



(0 



o (i) = 3 V w [lniV]. 



are the physical components of Eq.(|6|). This Euler equation is the version in physical compo- 
nents (with the extra dissipative terms J-^) of Eq. (2.29) of Ref. 1531 (see also Ref. ])55|). 



IV. RELATIVISTIC TRANSPORT THEORY 



We shall consider the phenomenon of relativistic transport of massive and massless particles 
(hereafter radiation) within a dense medium. We shall follow closely the formalism developed 
by Lindquist |5J|, but we will not repeat the details here. The radiated particles of the specie 
'R' will be treated classically as point particles except when interacting with the dense medium. 
The interactions and its quantum mechanical effects will be ultimately translated as emission 
rates and opacity functions. The particles will be characterized thus by a four- momentum 

where dX corresponds to an affine parameter (for massless particles) or to the proper-time per 
mass unit (for massive particles). 
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According to this formalism, one postulates the existence of a scalar function F-^{x^ ,p^) 
(the invariant distribution function for the specie R) which is a function from the phase space 
coordinates {x^,p^) to the reals. Actually, since we will be interested in particles satisfying 
the mass shell condition 

where m = 0, 1 for massless or massive particles respectively, will be a function defined on 
the reduced phase-space. 

The invariant distribution function so introduced will be such that the number density 
four- vector and the energy-momentum tensor of the particles are respectively [Eg] , 



j£ = J^F R (x^p^)dP . 
Tjf = J p^p v F K [x^jf)dP 



where dP is the invariant volume of the momenta space on shell [56]: 



dP = -NVh^ . 

Pt 

where d 3 p := dp 1 dp 2 dp 3 represents a coordinate 3- volume element and h is the determinant of 
the 3-metric. 

The distribution function F& is related to the dimensionless distribution function by 

Fr = ~ ^ 3 Fr ■ 

where <?r is the statistical weight of the particles of the specie R (e.g., g = 1,2 for neutrinos 
and photons respectively). 

As in the case of the perfect fluid, it is useful to refer the components of four-momenta to 
a tetrad. In the context of the 3+1 formalism we have, 

e : = p(*) = Np l , 
16 



with the inverse relationships given by 

P l = e/N , 
i N " i (i) 

where are the physical spatial components of the 4-momentum (i.e. the spatial physical 
components of the projection of p^ onto S^: 3 := h^p"). The ratio p^'/e corresponds to 
the local velocity of particles measured by Oe- 

Introducing the magnitude of the three-momentum as, 

2 (VI 

P -=P(i)P K1 , 



it is easy to see that (51) simply becomes 



2 2,-2 

e = p + m 



Here e is the energy (per mass-unit in the case of massive particles) as measured by the 
Eulerian observer. Therefore, from Eq. (|54|), one obtains, 

dP = Vh . d3p (l) , 

where we used that = ef^p^ = iVW) — pn\ef N" 1 = — iVe — pmef N l . 



When changing variables in the momentum space using (58) and ( pT|) a straightforward 
manipulations show that 

, 3 dp^dpWdpW ( pme^W 
dp = 7= 1 + A7 — 

where it was used the fact that detf^J = 1/yh. Then finally 

dp _ dpMdpWdp® 
e 

which has exactly the same form of flat spacetimes. 

The use of physical spherical variables in momentum space leads to 
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dP = p 2 dpdfl p /e = pded£l p . 

Indeed, this is the useful expression when dealing with spherical symmetry. 

The use of tetrad components allows us to write Eqs. ( |52| ) and (53), as follows, 

jgO = I p M FR{x \ p W) d p , 

T<f )(l/) = J p^p^F R (x\p^)dP . 
Therefore, the corresponding 3+1 matter variables are 

E R = J e 2 F R (x\p^)dP , 
J« = J ep^F R {x\p^)dP , 

= j ' p^p^F R {x\p^)dP , 
Sr = J P 2 F R (x\p^)dP . 
According to @ the total 3+1 matter variables are, 

E = £pp + E R , 

t(t) _ tW , tW 
— Jpp -t- , 

<?MC0 _ cWO') , oWO') 
j — °PF "i" R ' 

S = Spf + s R . 



where we remind that the quantities labeled with 'PF' are given by Eqs. (40)— ([43]). It is 
understood in these expressions that the sum of the different quantities extend to all the 
species R considered. 



A. The Boltzmann equation in curved space-times 

We defined a "macroscopic" 4-current density number of particles of the specie R by Eq. 
(p^). The number density of particles as measured in the local frame of an observer O p with 
four- velocity v^) is 

18 



rr 



R „, -n 



P 



J v^F R {x\p^)dP 



Therefore, the number of particles dN R with momenta between p^ and p 1 * + dp^ crossing the 
volume element dV of the space-like hypersurfaces orthogonal to and which is centered at 
some point x^ of spacetime is 

dN R = F R (x\pW)(-v^)dVdP . 

The quantity dW = {—v^p^) dVdX represents the four-volume spanned by the flow of particles 
(world lines) crossing dV, which is given by the element of hypersurface dV with normal 
and the particle's infinitesimal displacement orthogonal to dV given by dl = v^p^dX p6j . The 
quantity (—v^p^) dV is in fact the correct Lorentz invariant four-volume element. From the 
relativistic form of Liouville's Theorem (see Ref. |5(| for the details) dWdP remains invariant 
along a given set of trajectories. Therefore, the change in the number of world lines within 
dWdP is proportional to the change in F R 



5(dN R ) 



dx a + dy a 

8F R a dp a 8F R 
P + 



(-v^) dVdP 
dWdP . 



dx a dX dp Q 

The evolution for p^ will be thus governed by the equations of motion of individual particles: 



j\ ~~ P P 1 av ^ ^fields ^ •'coll 



This equations shows that the acceleration of the particles is due to 1) the space-time 
curvature 2) the forces arising from the interaction of the particles with fundamental fields 
other than the gravitational one, 3) the interaction with other particles that can be represented 
by "collisions". For our purposes, we will consider that J 7 ^^ = 0. That is to say, the only 
fundamental field we consider is the gravitational one. All other interactions like the weak ones 
(in the case of neutrinos) will be treated phenomenologically as collision terms, and therefore, 
the set of equations will not include gauge-fields, but rather involve macroscopic quantities 
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that characterize the medium and which are obtained from field theory in a similar fashion as 
one obtains the equation of state of the matter. 

The relativistic Boltzmann equation (RBE) then reads, 



LF R 

where 



coll 



is the relativistic Liouville operator often written by the fuzzy notation p a D/dx a (the direc- 
tional derivative of Fr along the phase flow), and (4dt) u = ~~ •^Rcoii^p" represents collec- 
tively the scattering, absorption and emission processes between the particles of the specie 
R and the medium. In the absence of collisions, the distribution function remain constant 
along the particle's path (i.e., along particle's geodesies). In the language of differential ge- 
ometry, the operator D/dx a = — P^aX'Spt corresponds to a coordinate basis vector of 



the horizontal part of the tangent space of the bundle B over the spacetime M 4 [56|. It is 
important to emphasize that the mass shell condition implies that the distribution function 
will not be defined over the entire tangent space but only on that part where p^p^ = —in 2 , 
that is the distribution function will be defined only on the sphere bundle (the subbundle of 
tangent vectors of fixed length). One can incorporate this restriction in the Liouville operator 



by treating the spatial part p l of the four-momenta as independent components. Then 



9FR_ v x T i 



B. Tetrad representation and the 3+1 RBE 



It will be convenient to use a tetrad components to re-write the RBE. By employing the 
tetrad formalism it's easy to show that the RBE reads 
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(oWn" — - „(0) V MO {S) 9 } Fx>(x a V^) ~ ( — \ 



coll 



Lp 



Let us consider Eq.(j33|), and split it in terms of temporal and spatial contributions. First 
we define, 

L p .-p p U {b){a)g ^ d) -p p ^ {b){a)de +V P U (b)(a) dp (i) > 

where the notation e = was used. The properties of the RRC and some straightforward 
calculations show a useful relationship between the four-RRC and the physical components of 
3+1 variables. For instance: 

u (tm ~ u (t)(t) ~ °M ' 

(t)(„) = °(0(m) = ( no sumation on m) 

U W(i)~~2l~ N 9(m) N wl + 9 (*rW e i W< ► UJI 

/T)(0 _ 3/^(0 

Here , are the 3-Ricci rotation coefficients, i.e., the RRC associated to the local basis 

WW ' ' 

frame on S t , and d {t) = £g + ^d {i) . 
In this way we obtain 

(«*<%•) - p^K {m ) §- + (e 2 a» + _ + p (0 p W 8^) jL 

Finally, the 3+1 decomposition of Eq.([33"l) is, 



e%) + - (ep^a {i) - p (l) p^ K {l){j) ) 



coll 



It is to be noted that the properties of RRC imply that ^('vn (no sum convention) are null. 
Therefore, the above equation can be further simplified. 

Equation (|87| ) is the 3+1 version of the of RBE, here written in terms of physical compo- 
nents. The mass shell condition e 2 = p^pu) +fh 2 can be imposed on the RBE by considering, 
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for instance, pW as independent variables. In that case is to be considered as though it does 
not depend explicitly on e, i.e., dF^/de = 0. Alternatively, the use of spherical-like variables 
in momentum space (see next section) will allow us to consider e as independent variable and 
p 2 = as the dependent one. 

C. Collisions 

As we stressed before, particles may be submitted to collision forces arising from the inter- 
acting medium (in the case of neutrinos these forces come from weak interactions with baryons). 
Let us remind that the collision integral as conceived originally by Boltzmann assumes that 
the interacting medium has a known distribution function. That's it, the distribution function 
of the medium is a data of the problem. 

In the present case, the interacting medium is to be considered not as particles but rather 
as a fluid field, namely a perfect fluid. Thus, for our purpose, it will be more convenient to 
characterize the collision integral in terms of scalar functions as it is usual in transport theory. 
In this way collisions will be represented macroscopically by the so-called invariant opacity 
o(x lx ,p tJ ') and the invariant emissivity T(x fl ,p fl ). 

We then assume that the collision integral takes the following form 

( d\ ) co u R 

In terms of quantities measured in the same frame one can write 

T _ n{x^p») _ n'jx'^p'n 
e 2 e' 2 

here r\ and e being the matter emissivity and the particle energy respectively, both measured 
in the same frame. In the same way, one can introduce the matter opacity as 

o = e X (x^pn = e'x'(x'^p'n , 
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where x ~ 1/^, ^, being the mean free path of the particle in the corresponding frame. 
The collision term takes then the useful form 



(dF R \ 
\ dX / coU 



XF 



The opacity x an d the absorption coefficient k are related by 



X = KU , 

where n is the proper number density of particles that composes the medium (e.g., the baryon 
density), such that n := —j^u^ where is the four-current of baryons. 
In this way, an alternative form of the collision integral is 

(^L = ^ ( *- Fr) ' 

where k e = en, and S = T/o is usually referred to as the effective source function. 

It is to be emphasized that quantities measured in different frames are related to each 
other via the invariant quantities and Lorentz transformations, for instance, the relationship 
between the opacities meaured in the Eulerian frame and those of the proper frame of the fluid 
are given, according to Eq. (|90|) , by 

ex = dpXp , 

where quantities in the left-hand-side (l.h.s) refer to the Eulerian frame, while the quantities 
in the right-hand-side (r.h.s) refer to the proper frame of the fluid. Since physical components 
of four vectors in both frames are related by a Lorentz transformation, for instance 



where 



( A M P 

with 3 [/W given by Eq.(| 



V 



p. 



r - 3 £/»r 



and r by Eq.(44), then for the time components, 
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e p = r (e - 3 U (i) p {i) ) = eT (l - 3 tf«t; (i) ) = eT (l - || 3 C/ (i) || cos(# R )) , 

where t)W := p« /e represents the velocity of the particles with respect to the Eulerian frame. 
In the last formula one recognizes the well known formula for the energy shift due to the relative 
motion of observers. The type of shift (red or blue) will depend on the angle #r between the 
propagation vector of the particles j)« and the velocity of the fluid 3 LA*' (i.e., blue or red shift 
if the fluid is approaching or receeding respectively from the Eulerian observer). Therefore the 
transformation formula between opacities yields 

x = x P r(i-\\ 3 u^\\\\v^\\co S (0 R )) . 

In the case of massless particles | \v^' \ | = 1. In a similar way one can obtain the transformation 
formulae for the absortion coefficients and the emissivities. 

D. Conservation equations for the radiated flow 

The particle number current and the energy-momentum tensor of the radiated particles 
was introduced by Eqs. (|52|) and (|53|) respectively. For instance, in the case of perfect quantum 
gases in thermal equilibrium (i.e., Fermi and Bose gases) the above definitions allows one to 
recover the usual macroscopic expressions for the energy-density, density-number and pressure 
parameterized by the temperature, particle-mass and chemical potential of the species as 
measured in the local frame. 

When collisions are present, both the particle number and the energy-momentum tensor 
of the particles will not conserve alone since there will be exchange of energy and momentum 
with the interacting dense fluid. Thus we can expect that the conservation equations derived 
from Eqs. ( |52|) and ( |53|) will have sources arising from the collision integral: 

and 
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Then we write 



where 



V„j R = K R , (101) 
V„Tg = -J* = nv& , (102) 

K K := J n(x^jf) [S{x^jf) - F K (x^p^\ dP , (103) 
w ^ : = J p» K {x»,p») [S(x",pi") - F R (x^ l ,p' 1 )} dP . (104) 

One can define the mean emissivity (energy/ (volume x time) in the fluid frame as 

V := -u^V^ = -nu^ . (105) 

Since we have been using quantities measured in the Eulerian frame, in the above expressions 
we are to use the corresponding quantities with respect to the same observer. Namely, the par- 
ticle number density tie as measured by the Eulerian observer is related to the proper number 
density of the perfect fluid n by tie = nr. Same considerations apply for the remaining collision 
variables. For instance, the emissivity measured in the Eulerian frame T>e = —n^nw^ = nw^ 
is related to the proper emissivity V = —u^nw^ by V = TT>e (l — 3 U^ where 
:= 



V. THERMODYNAMICS 

In this section the thermodynamical description of the dense matter with which the ra- 
diative particles interact will be presented taking into account the mean quantities introduced 
in the previous section. Such a description is performed in the proper frame of the fluid. We 
then assume that the equation of state of the dense matter (i.e., the perfect fluid) is given in 
parametrized form as follows 
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p = p(s,m,... ,n m ) , 
p = p(s,ni,...,n m ) , 

where s is the entropy density and nu (1 < M < m) is the number density of particles of the 
specie M (e.g. baryons and the different lepton flavors), all of them measured in the fluid frame. 
For instance, in the case of a dense matter in hydrostatic equilibrium composed by a mixture 
of neutrons, protons and electrons, the electron density is obtained directly from the proper 
baryon density n by demanding charge neutrality and chemical equilibrium p n = p p + p e . This 
last condition arising from the equilibirum of the nuclear reactions: n ^ p + e~ [Sj. In that 
case the equation of state depends on n solely. 

Now, the equations ( |106| ) and (107) are not independent from each other but linked through 



the first principle of thermodynamics 

dU = 9dS - P dV + / u R dN R 

where and ^ R are the temperature and the chemical potential of the specie R of the particles 
composing the fluid respectively, defined by 



6 



ds) 



» R 



Using these definitions, Eq. ( |108| ) takes the following form in terms of densitized quantities 

p = 0s + /u R n R - p . 

This equation is often referred to as the compatibility themodynamic condition between Eqs. 



pa) and (107) M. 



The conservation equation for the baryon number reads, 

V^ = . 
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where we remind that = nu^ is the density current of baryons and n the proper baryon 
number density. 

This equation can be written explicitly as an evolution equation for the number density 
he '■= —n^j^ = nT measured by the Eulerian observer as follows 

d t (Vhn E ) + di IVhriEV 1 ] = , (113) 



where T and V' L are given by Eqs. ( f44| ) and ( f46| ) respectively. Introducing the physical 
components of the fluid velocity field given by Eq. (|4^), we have the alternative expression, 

3 t (7W) + di [Vhn E (N* + Nqfa 3 U^)] = . (114) 

The equation of conservation of baryon number leads to the conserved total baryon number 
given by 

M= [ -fn fl Vhdx 1 dx 2 dx 3 (115) 

= { nTVhdx 1 dx 2 dx 3 . (116) 
JT, t 

The integral has compact support corresponding to the volume enveloped by the star surface. 
In a similar way, the equation for entropy conservation reads 

V„0m") = ~(M B fcR + Z>) , (117) 

where 

K R := V M (n KU i*) , (118) 

is the rate of particle production and we remind that 

V := -<VvVC = u„V M T^ , (119) 

is the particle's mean emissivity in the fluid frame. Therefore the sources for the entropy 
generation in a perfect fluid are from the particle production (e.g. neutrinos). The Eq. ( |118| ) 
is completely equivalent to Eq. ( |101|) which is given in terms of the distribution function. 
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One can define the entropy per baryon a = s/n and use Eqs. (|112|) and (|117|) to obtain 



Explicitly this provides an evolution equation for a: 

d t a + V%a = -^-L R KR + v) . 

refc)l \ ' 

Moreover, using Eqs. (f44|), (^5j) and the tetrad approach of Sec. IIA, the Eq. ( |120|) takes 
the alternative form: 



nre 

In the same way, we can introduce the particle number per baryon = n^/n and write (111 



as 



u^V^xr = -IZr , 



n 



which provides an evolution equation for xr: 



or alternatively, 



N 

d t XR + V l diXR = —TZr , 
ni 



n^XR + 3 U^d {l)X R = -^TZr 



nT 

In order to close the whole system of equations presented so far one needs the input 
of particle physics. That is, the equation of state for nuclear matter, the rate of particle 
production and the opacities (see Ref. [|| for a review). In the case of neutrinos emitted by 
nuclear matter out of beta equilibrium via direct and inverse (3 processes during neutron star 
collapse, the rate of particle production, the emissivities and the opacities can be given in 



terms of rather simple formulae pq,p7| (see also Ref. [58| for neutrino emissivities from quark 



matter in (3— equilibrium within neutron stars and Ref. |56] for neutrino emission from hot 
and dense atmospheres). For the case of reaction rates and opacities in Type II supernovae 
see for instance Ref. |I5| . 
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VI. SPHERICAL SYMMETRY 



In this and the following sections we will focus on spherically symmetric spacetimes. The 
most general line-element for such spacetimes according with the 3+1 decomposition of the 
metric Eq. (||), writes 



-(N 2 - N r N r )dt 2 - 2N r dtdr + A 2 dr 2 + B 2 (r 2 d9 2 + r 2 sin 2 



(126) 



where all the metric potentials are functions of the coordinates r and t solely. The three-metric 



hij is easily read-off from ( 126 ). 

On the other hand, the triad coefficients are 



ey = diag (A(r, t), r B(r, t), r sm9B(r, t)) . 



(127) 



The inverse coefficients qhs can be obtained trivially form ( p.27| ). 



The extrinsic curvature can be computed from (|J). We find [30 
(d t A + S^JVW) 



/-A 









rB 2 / N (r) r _ 

"TV" {— + B d tB 




+ 



r d r B 
AB 



) 








rB 2 sin 2 (6>) 
N 



/-h( 



(*',) - (*%) 



d t A , N r d r A , 
A A 








d r N r ) 



N \ B 



i N r d r B , iV^ 



/ _J_ (dtA , ^JVW 



V 












1 fgtg , 
N \ B 





{— +B d tB+ 35 J / 



(128) 






N r d r B | TV 1 / 



B 



rA 







(129) 



J_ /dtB i jyW d r B i ATW \ j 
JV \ B ' AB ' rA J ' 



where the index of Eq. ( |129| ) was raised with h lJ from Eq. (|12 
The three-scalar of curvature is given by 



A 2 



A 2 f 2d r A 2{d r B)(d r A) (d r B) 2 6d r B 2d 2 r B^ 



A 2 \ rA 



BA 



B 2 



rB 



B 



(130) 



29 



A. The 3+1 Einstein equations 



It is useful to introduce the new variables, 



v := ln[N] , 

a := ln[A] , 

/3:=ln[B] . 

The Hamiltonian constraint Eq.(|i~0|) reads 



1 [£. 

r 2 Is 2 



1 +A 2 



(6) \2 



2d r a 6d r P 



+ 2(d r a) (d r p) - 3(d r p) 2 - 2 <9 r 2 r ./? = 8ttG EA 2 . 



where we have used the fact that the non-diagonal components of K l , are null and that 



K 



(9) 
(9) 



K% [cf. Eq. (HI)]. 



The radial component of the momentum constraints Eq. (|ll|) , reads 



,(r) 



K 



{ K (r) - 



+ d r (3 ) - d r K (6) {e) = 4 vrGo A J (r) , 



Or in terms of the trace K, 



K 



(9) 
(9) 



-+d r /3j -d r K + d r K\ r) 



8irG AJ, 



(r) 



The angular components of Eq.(ll) and the spherical symmetry lead to the conditions 
Jg = = Jfh which implies the absence of "angular currents". 



The dynamical equations for the non-diagonal components of Eq. (|i~2|) , i.e., for dtK^ r ' 



(6)> 



dtK {r ^ } , d t K {6) {<t)) with the fact that K (r L 



respectively to the conditions that S 



(r) 
(9) 



s 



K 



(r) 



(r) 



K 



(9) 

(4>) 



[cf. Eq. (129) ], leads 



S 



(9) 



0. Moreover, taking into account 



the fact that K^X\ = the dynamical equations for dtK^Xs and dtK^Xs [see Eq.(|T^)], 

leads to the condition S X\ = S^X\ corresponding to an "isotropic" energy-momentum tensor 
which is compatible with the hypothesis of a spacetime with spherical symmetry. In this way, 
the only two non-trivial dynamical equations are 
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iV« d r K (r ) 



" NKA _IL(3^_^ll + 2{dra) {dr(3) _ 2{dr(3 f _ 29^ 



° tn (r) + A ( r ) A 2 V r 

2V 



£> r > + (d r u) 2 - (d r a) (d r v)] = 4vrG N (-S {r) {r) + 2S {e) {e) - E) , 



d t K% + 



d r K {e ), 



A 



w NKK m 



N 

(<?) ~ 



+ + {dra) {dr(3) " 2{dr(3)2 " 



+ 



AT 
12 



(d r /3) (CU) + 



d r v 



AttG N (S ir) (r) -E 



From Eq. ( |l4|) the evolution equation for the trace of is 



d,K + iV( ' > 8rK - N 



A 



{K {r) (r) ? + 2{K^ 



(0) i 



+ 



N 
A? 



d 2 „v + [d T vf - [dra) (d r v) + 2(d r f3) (d r u) + 



2d r 



4irG N (S + E) 



where we recognize the 3-covariant Laplace operator in spherical coordinates of the slices 



S t [sec Eq. ([1261) 1: 



A 2 



8Lv + 



2d r 



(d r a) {d r v) + 2{d r p) (d r u) 



B. Matter equations 



In the present case of a perfect fluid in spherical symmetry we have 



where 



Epf = (p + p)T -p , 



PF S% = (E PF +p)( 3 U^) 2 +p , 



(r) 

<j(«) _ <?(*) _ _ 
PF J ,vn - o ia^-P 



(9) ~ » (fl 



'(r) 

r 



iV v ; N V 
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with V r := u r /u t . The spherical symmetry implies u 9 = = and therefore = = ?JW. 
Note that 

PFS%= 3 U^J^+p . (147) 

The Eq. © reads 

d t £ PF + iV^PF + (abVJpf) = N (pF^iT^ + 2 PF S (e ( ) e) i^ ( ) e) + £ pf k) - 2J£ p P iV - iV 2 ^ . 

(148) 



Or in terms of Ji r \ = AJ r , 

d t E PF + A^ r £ PF + (BV J pf ) = iV ( PF^if % + 2 PF S { d) {e) K { % } + E PF K) - | J$d r N - N 2 P . 



(149) 



Using Eqs. ( |143|) , (|146| ) and ( |147| ) the latter equation can be written in conservative form as 
d t E PF + \d r (r 2 V r E PF ) = E PF (nK + — + d r N r } - \d r (r 2 ^ 3 U^ p 



NJ, 



PF 



(r) 



A 



d r u + d r a + 2d r (3 - 3 C/ (r) AK (r ) r) ] + NKp - N 2 T l . (150) 
The momentum conservation Eq.([l9|) reads 
d t J r PF + N r d r J r PF + J? F d r N r + N{ 3 V ; PF <S' r ) = NK J r PF - ( PFS r r + £ PF ) d r N - 3 F r N . (151) 
Explicitly 



dt J PF + N r d r J r PF + J PF <9 r iV- + Nd r PF S {r l + 2^V f p F S ( l - pfS { 1 ) I - + 



1 d r B 



= - ( pfS^ + £ PF ) 3 r iV + AT J r PF (i^ + 2 if (( jy - 3 ^ r iV . (152) 
Or in terms of Jm, we have 



'W ^ » "X ^ " r "W ^ "F 1 ^ " ^ T Ur PF " W ^ "X V PF " (r) - ^ wj y- + "b 
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pf^W , + £ PF ) M + at jpf (jf W ) + 2 jf W ) _ 3^ )iV (153) 



When using Eq. ( |129| ) to replace the time derivative of A we find 

8 rPF , ^ (r) a jPF.N (r) 2N , (r) (e) x /l d r B\ 



F {r) N . 



Using Eqs. ( |146| ) and (147), one obtains an equation for jff in conservative form 



AT 



PF 



3 C/ (r) {d r a + 2d r (3) - 4 (~N r + drN 1 
N \r 



(E PF + p) d r v - d rP - 3 F(r)A\ 



with the alternative form for the r.h.s, 



d t Jf r f + \d r (r 2 V r J$) = J ( p r f \2N (K {r l } + *T%) - V r {d r a + 2d r f3) + AT ( <7 ; n + 2^. J - - ) - rV,..V 



iV r 



(Ep F + p) d r u + d r p + 3 J r {r) A 



Another possibility which has turned to be very useful in some numerical studies P|9|,|36|] , is 
the use of the Euler equation for the fluid instead of the equation for J^ r y The only non-trivial 



component of the Euler equation (45) in spherical symmetry reads 

1 



d it ) 3 U^ + 3 U^ 3 d M 3 f/W 



(r) 



E FF +p 

+ ^(3^. Al » ) _ aw)+ _i_ ( . Wl ,_^, ) . 

(r) (r) 

where we used the fact that the RRC are antisymmetric so that the terms ®u)M = ^ = (r)(r) 
and we remind the explicit expressions 



1 » JL N ''' 

1 



<%) = — d t + — d r , 



d( r ) = -jd r , 



«(r) = 3 d^U 
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C. The RBE in spherical symmetry 



The most general 4- metric for a spherically symmetric spacetime is given by Eq.( |126| ). 
Then, it's easy to see that the simplest tetrad choice associated to Eq.( |126j ) and which 
corresponds to the local tetrad of the Eulerian observer, reads, 

e w = iv^ + lv^' (161) 

e « = il' (162) 

e W = VBd9' (163) 

e ^ = 7B^6d4> ■ (164) 

where d/dx^ [with x^ = (t,r,6,(j))] denote the coordinate basis. 

In the Eulerain frame, the spherical symmetry of configuration space (i.e., spacetime) 
induces a symmetry on momentum space that can be exploited to simplify the computations. 
It is convenient to define spherical variables on momentum space as follows: take a unit vector 
e( r ) as polar axis, i.e., as symmetry axis on momentum space. Then it is useful to introduce 
new variables e, 0, 7 in the Eulerian frame as 

p® =e, p^ =pcosifj, p^' = p sin -0 cos 7 , p^ = p sin -0 sin 7 . (165) 

As mentioned earlier in Sec. IV, e is the energy of the radiated particles (e.g., neutrinos) in 
the Eulerian frame and p 2 = p^'pu) ; ip is the angle between the polar axis and the neutrino 
propagation three- vector pW, and 7 is the angle of rotation around e^ r y The above variables 
are consitent with the mass shell condition (|6T|). For massless particles, then p = e. Under 
these new variables the RBE (|83l) reads 



coll 



where we have explicitly stress the dependence of momenta with respect to the new momenta 
spherical-like variables (e,p,tp,j) represented colectively by p. 
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(166) 



Indeed, the spherical symmetry and the mass shell condition will be reflected in the RBE 
by the fact that the distribution function F R depends only of four phase-space coordinates 
(i, r, e, //) or (t,r,p,/j,) instead of the original eight {x^^p^) |]56| , |6l| , where := cos^ • 

We can now proceed to calculate explicitly the RBE. The only non-null Ricci coefficients 
are [|6C|1, 
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U WW " U toW~ rB ' 



where := -^Jj , D r := -^J^ = 3 <9( r ). After imposing the mass shell condition, the RBE 



(|166D in terms of the spherical variables reads explicilty 



(pa A/"W\ 
eD t F R (r, t,e,fi) + l — + 1 eD r F R (r, t, e, /i) 



(1 - ^ 2 )^ {— A + D r p) + M 2 Aa + (1 - li 2 )D t f3 + ^D r v + ^-D r N^ 
I\ \rA J p TV 



+e(l - /x 



/ co ii 
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N I \rA J p JS 



d e F K {r,t,e,n) 
d^F K (r,t,e,/j,) 



The alternative 3+1 form of the RBE can be computed from fl87| ) when changing to the 
spherical variables in the momentum space: 



(vu iVW\ 

eD t F R (r, t,E,n) + l^ + —\ eD r F R (r, t, e, n) - p 2 



(r) 



P 



d e F R (r,t, 



e 

—D r u 
P 



dF\ 



In normal coordinates where N r = 0, and for massless particles, Eq.( |166D reduces to the 
relativistic Boltzman equation derived by Lindquist E3| with the choice B = R{r)/r. Under 



the isotropic gauge Eq.(168) can be written in a "conservative" form which is specially suited 
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for numerical solutions [45|. Following Harleston & Vishniac [44], we can write Eq.(168) in 
conservative form as 



o t F + ^d r (t 2 v;f) + jd e (fn e F) + ^ (h^f 



NAB 2 fdF\ 
~d\) 



coll 



where 



F = AB 2 F R (r,t,e, f i) , 

! „ , ,1 + N 



p l A \ e 
W e = -L(e) = -iV 



-L(p) = N(l - v 2 



.(l-f)Kl%-fK% + !fD r v 



e 

—D T v 



e \rA J \ v i vis p 

and L stands for the Liouville operator as it appears in the l.h.s of Eq.( |168| ). Alternatively, 
the above equation can be written as 

0,F + ^ (rXP) 4- i 8p ( P °H r F) + ffc («„f) = ^ , 

where now F = F&(r,t,p, fj) is to be regarded as a function of p instead of e and 



N - 

n p = —L{p) = -n 

pe 



2^(0) 



(0) 



When one of the different gauges and slicing conditions are chosen to write the RBE this 
will only affect the particular form of the quantities given by ( |170D — ( |173| ) and ( |175| ) , and the 



r.h.s of Eqs.(|i69Q and (g7g). 

Finally, we emphasize that all the momentum variables which appear in the various forms 
of the RBE in spherical symmetry, are components with respect to the orthonormal tetrad 
carried by the Eulerian observer. Therefore, the corresponding quantities in the collision 
integral are to be referred to the same observer. 
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D. Mean radiative variables 



In Sec. IV we defined the energy-momentum tensor of particles in terms of their microscopic 
four-momenta. However, we can introduce mean radiative variables which are to be interpreted 
as their counterparts of the continuum case. Such variables are called the moments of the 
distribution function. In order to write them explicitly, we shall use the physical components 
of the energy-momentum tensor of particles Eq. (^) and the spherical variables in momentum 
space in addition to the invariant volume element of momentum space given by Eq.(^). We 
have then 

coo r\ /•27r 



t (m)W = f°° f f n p Mp(v)F n (r,t,p,u)—dpdud7 
Jo J -i Jo e 



OO rl /•2-7T 

p WpWF R (r, t, e, n) \Je 2 - rh 2 dedfidj . (176) 



fh J-l JO 



E R := TjP = 2vr / jfyjp 



The only non-null moments are 

/•oo i r-1 poo rl 

/ p 2 \/p 2 + rh 2 dp / Fn(r,t,p, ii)d/j, = 2ir / e 2 Ve 2 — m 2 de / F R (r,t, e, n)dfi , (177) 

Jo J-l Jfh J-l 

H R := Tjp } = 2tt r P 3 dp f 1 fxF R (r,t,e,fx)dn = 2vr fe(e 2 - rh 2 )de C fiF R (r,t,e, n)dfi , (178) 

JO J-l Jfh J-l 

p R := 4 r)(r) = 2vr P dp f 1 p 2 F R {r,t,e,fi)dfi = 27r f°° (e 2 - m 2 ) V2 de [* fi 2 F R {r,t,e, fi)dfi . (179) 

JO yjp z +m z J-l Jfh K ' J-l 

The tangential preassures are given by, 

pi := if W = Tf W = vr f°° -j^—.dp f (1 - » 2 )F R (r, t, e, »)dp 

Jo \/p A + m z J-l 

f°° / \ 3/2 

= irj [e 2 -m 2 j de J (1 - fi 2 )F R (r, t, e, {i)dp (180) 
We note that in the massless case, 

pl = \(E R - PR ) . (181) 
The corresponding 3+1 variables are 
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J R r) =H R , (182) 

s£ )(r) =PR, (183) 

S^=S^=pl , (184) 

Sr = r,S% = Pr + 2p£ . (185) 

The effective pressure of radiation which can be defined as p^ s = Sr/3 turns to be in the case 
of massless particles p^ = E R /3 which corresponds precisely to the equation of state (EOS) 
of an ultrarelativistic gas. In Eq. ( |178|) , Hr is the mean radiative flux of energy in the radial 
direction. 

It is usual to introduce the so called variable Eddington factor 



Jo J-iJo e 



Ifh J-UO 

In particular, the mean number density and flux of particles measured by the Eulerian observer 
are given respectively by 



ra| := -?V7r = 3r = 2vr J J ^ F R (r, t,p, fx)p 2 dpdfi 



(186) 



„ _ PR 
^R 

used to measure the degree of "anisotropy" in the particle flow. In the case of massless particles 
if S = 1/3 then p^ = p R = Er/3 which corresponds to an fully isotropic flow. Moreover, in 
the free streaming approximation we have 

Er = Jr T) = 4 )W , (187) 

and so S = 1, which is the case of a highly anisotropic flow (p^ = 0) with a purely radial flux 
of radiation. 

In the same way, the macroscopic particle number density current measured in the Eulerian 
frame is given by Eq. ( p6| ) and it in terms of the spherical variables of momentum space gives 

coo rl r2iv . . n 2 



oo 1*1 p2n 

p M F R (r,t,e,n)Ve 2 -m 2 dedfid-f . (188) 
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/■OO /*1 /•27T 

= 2ir / / / eFR,(r, t, e, / u)Ve 2 — m 2 dedfi , 

im V-l JO 

3r = 271 i o - o <fo / fJ>F R (r,t,e,n)dfi = 2Tr (e 2 - m 2 )de fj,F R (r,t,e, /j,)d/i , 
Jo \/p z + m z J-i Jfh J-i 

andjW=0 = jf . 

In terms of the above macroscopic variables, the energy- momentum conservation equation 
in spherical symmetry Eq. (|100| ) reads according to Eqs. (36) and (37) as follows [60|: 



d (t) E R + 3 8 {r) H R - E R K + 2H R 



d(t)H R + 3 d< r) p R + 



^4 + ^O/ 3 



[3m - E R ] + [p R + E R ] 3 d {r) v -H R (K + K^) = -F® 



Depeding on the gauge and slicing choice some of the terms within these equations can 
vanish. Actually such evolution equations can be obtained directly from the RBE when mul- 
tiplying this by the momenta and then integrating in momentum space. 

We emphasize that it is more convenient to calculate E R , H r and p R , directly from their 
definition once the distribution function has been computed, rather than using the above 
equations. In any case, such conservation equations can be used to verify the self-consistence 



of the system. The disadvantage of using the system of Eqs. ( |191[ ) and (192) for the moments 
of the distribution instead of solving the RBE is that such a system is undetermined (i.e, there 
are more variables than equations). Then a closure relation is needed to remove the ambiguity 
(e.g. the diffusion approximation relating E R and p R ). 



(189) 
(190) 



(191) 
(192) 



d t (AB 2 n E ) + \d r (r 2 V r n E AB 2 ^) = 



d t (AB 2 n E ) + \d r 



AB 2 r 2 n E (N r + ^ 3 U^ 



A 



ner 
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Finally, in spherical symmetry the evolution equations ( p.!3| ), ( |114| ), ( |121| ), ( |122j ), ( |124| ), 
(|125D write respectively 



(193) 
(194) 
(195) 



N 3 U^ N 

d t a + N r d r a H d r a = — 

A nL<8 



^TZ R + V p \ , 



N 

d t XR + V r d r x R = —TZr , 
nl 



d t x R + N r d r x R H d r x R = —TZ R . 

A nl 

Equation ( |193| ) has a conservative form for the quantity AB 2 ue- One can alternatively 
write an equation in conservative form for n E as 

d t n E + \d T (r 2 V r n E ) + n E [d t a + 28 t (3 + V r (d r a + 23 r /3)] = . 



In a similar way, 



OtnE + 4r<9, 



N 



+ n E 



8 t a + 28 t p + ( N r + ^ 3 U^ ) (d r a + 2d r /3) 



. 



The total contribution of sources [Eqs. d72|)— (|75[)] that appear in the 3+1 Einstein equa- 
tions write in spherical symmetry as follows: 



E = E R Y + E R , 

j(r) = AJ r = (jBpp +p) jj{r) + Hr ^ 

S (r l ) = (E PF +p)(U^) 2 +p + p R , 
s (9) _ s (4>) 



(9) - w (<f>) - P + PR ' 

w 



S = S { % = (E PF + p) {U^f + 3p + PR + 2p£ . 



VII. GAUGE AND SLICING CONDITION 



Two of the most popular gauges and time slicings that have been used in spherical symme- 



try are the radial and isotropic gauges and the maximal and polar slicings (see Refs. [ ^9|j62 ,63] 
for a more general discussion on the slicing choices). Moreover, it is in the framework of asymp- 
totically flat spacetimes (condition usually demanded in astrophysical applications) that they 
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become specially useful; it is in this context that such gauges and time slicings will be discussed 
in the following. 

The isotropic gauge A = B with the maximal slicing condition K = 0, dfK = (IGMS) 



case, that choice leads to the well known Schwarzschild solution in isotropic coordinates. The 
maximal slicing has the advantage of freezing the evolution in regions near the formation of 
space-like singularities while allowing a faster evolution in the outer regions (feature usually 
quoted as "singularity avoidance" property). The time slicing leads to an elliptic equation for 
the lapse and therefore, for rather general matter conditions (e.g. strong energy condition), 
one can use the maximum-minimum principle to determine the qualitative behavior for the 
lapse (cf. p9| , |62| ). The lapse function has a minimum at r = and a maximum at r — ► oo and 
during the evolution the minimum tends to zero as t —* oo (the collapse of the lapse) halting 
the propertime separation between neighbouring slices as the singularity forms. However, 
far from the origin TV — ► 1 which allows to advance the evolution in the asymptotic regions. 
The IGMS coordinates have also the advantage that the three-metric remains regular at the 
formation of apparent horizons, which allows to continue the evolution. The drawback is that 
eventually the metric potencial A grows exponentially at the origin and then the coordinates 
are "sucked down" to the black hole, which avoids a good description of the evolution outside 
the event horizon. 

The radial gauge B = 1 with the polar slicing condition (RGPS) K = K r r has been em- 
ployed sistematically by Gourgoulhon [||,||J36| . These coordinates are a generalization of the 
Schwarszschild coordinates to the non-static and non-vacuum spacetimes. The field equations 
turns to be much more simpler than those of the IGMS coordinates since the equations for A 
and N reduce to first order in r while evolving in time through their sources. Furthermore, 
the shift N r is zero everywhere on the slices. The RGPS coordinates has a central "singu- 
larity avoidance" property which is even stronger that the IGMS coordinates. In fact, the 
slowing of the evolution is such that it avoids the formation of apparent horizons, the metric 



has been employed by several authors 
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potential A, however, diverges at the star's surface as the matter enters the Schwarzschild 
radius §§,|lj]], leading, unlike the IGMS to a coordinate crash. Thus, these coordinates do not 
serve to describe the black hole interior. Nevertheless, the pathological behavior of the RGPS 
coordinates occurs at a large t, and from the astrophysical point of view (e.g., from the point 
of view of an observer at spatial infinity), these coordinates are good enough to describe the 
entire evolution of matter outside the black hole, the ingoing matter takes by the way, for an 
observer at infinity, an infinite time to cross the event horizon (the evolution is thus "frozen"). 

The hybrid choice of isotropic gauge and polar slicing (IGPS) has been less used in the 
past. However, it seems that they overcome the drawbacks of the above coordinate choices (cf. 
1 1 , 4C ] ; see also Ref. [63] for an analysis of these coordinates in the context of axisymmetry) . 



Another popular gauge choice is the comovil gauge (Lagrangian gauge). This gauge has 



been particularly used in the study of supernovae collapse with syncronous [13,24,E>| pH] and 



polar slicings |7jBT|. The comovil coordinates have the advantage that the hydrodynamic 
equations are simpler since there is no advection. The synchronous slicing are orthogonal in 
the sense that there is no shift (N r = 0). The disadvantage of the latter is that they fail badly 
when black holes start forming (cf. 0,^]). The asynchronous coordinates can remedy this 
problem. In particular Schinder and coauthors [f7,|61|, have used polar slincings to avoid the 



pathologies of the synchronous gauge. Another modification of the comovil and syncronous 
coordinates that allows to handle the formation of black holes is the introduction of an outgoing 
null coordinate instead of the usual time coordinate f|27|J6£ 



Finally, we mention the isotropic gauge and constant-mean-curvature slicings 
(IGKt— slicings) K = K(t), employed by Harleston and coauthoros [ 44 , 45 1 . This time slic- 



ing contains the maximal slicig as a particular case. It also posseses the feature of strong 
crushing coordinate avoidance. These coordinates generalize (to the non-homogenoues case) 
the comoving coordinates of homogeneous and isotropic spacetimes which are relevant in the 
standard cosmology. Such a choice is thus useful when the space-time is required to be asymp- 
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totically Friedmann- Robertson- Walker [45,01]. The main difference between this choice and 
the maximal slicing condition is thus the behavior of the hipersurfaces asymptotically: the 
maximal hypersurfaces reach spatial infinity while the -fQ-hypersurfaces reach future or past 



null infinity whether K is positive or negative [37 



A. Isotropic gauge and maximal slicing condition (IGMS) 



The isotropic choice A = B (a = (3), implies from Eq.(|129|) that 



K r „ = K {r) 



(r) _/V 



(dta + d r N r + N r d r a) . 



Therefore 



3 t a = —NK r - 8 r N r - N r d r a . 



On the other hand, the maximal slicing condition K = implies that 



K\ = -2K\ , 



or equivalently 



3d t a + 3N r d r a + d r N r + 



2N r 



Using (207) in previous Eq. we obtain 



d r N r 



-NK r 



This can be written as to give the following differential equation for N r [cf. Eq.(21) of Shapiro 
& Teukolsky (1980) §], 



Or 



N r 



2r NK 'r 



On the other hand, using Eq. ( 210 ) in Eq. fl207|) one obtains an evolution equation for a: 

d t a 
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-n* Ua + 1) + Ink* 



With the above choice and with fl208|), the Hamiltonian constraint Eq. (|134|) reads 



2d z rr a + {d r aY + 



-8ttG EA 2 - jA 2 (K r r ) 2 



Adopting the variable 



a = a/2 , 



we obtain a second order differential equation for d 

2 ad 



d 2 r a + 



2nG E + ^-(K r r ) 2 
lb 



(ad) 2 , 



where we recognize in the l.h.s the Laplacian operator of a spherically symmetric Euclidean 
space [cf. Eq.(19) of Shapiro & Teukolsky (1980) || for a source term E including a perfect 
fluid alone Eq.flMl])]. 

The momentum constraint Eq.(136), reads 



3K\ \- + d r a)+ d r K\ = 8vrG J r , 



were we used J r = AJt r y This can be written as a differental equation for K r r as [cf. Eq.(20) 
of Shapiro & Teukolsky (1980) §§], 



a U 3 r 3 K r 



\irGnr 3 A 3 J r . 



The Eq.( |139|) provides an elliptic equation for N [cf. Eq. (18) of Shapiro & Teukolsky 
(1980) § for source terms E and S of a perfect fluid Eqs. ( |i4l|) and Q], 



,2„ , ,,,,,2 , ,;,..,,,■;,„, , 2 ^ = 47tGqA 2 {S + E) + 3 A 2 {K r )2 

r 2 



9 2 r ^ + (d r u) 2 + (d r a) {d r u) + 



We have then four differential Eqs. (|T|), (|T|), (|216|) and (|218|) for iV r , ,4 and ET r r y iV 
respectively. It is to note in those equations that the field variables evolve in time through 
the matter fields. Although the evolution equation for K r r is redundant, for completness we 
write it in the IGMS coordinates. Equations ( |137| ), ( |213 ) and ( |218j ), lead to 







+ A + 4 I W) A' 



(d r u + d r a) ( — I- 9 r a ) + (d r u)(d r a) 



SttGqNS 



(r) 
(r) 
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Concerning the energy conservation equation ( |148| ) , this reads, 
d t E PF + N r d r E PF = -Nd r J PF - NJ PF \ 2d r v + * + 3d r aj + NK {r) {r) ( PF 5 (r) 



(r) 



PF S W {0) ) - N 2 ^ 



Or in terms of J r = A 2 J r we obt 



am 



d t E PF + N r d r E PF = -^d r J^ - J r PF [2d r v + -+ d r a) + NK^ ( PF S^ r) - PF S^ e) ) - N 2 ^ 



The conservative form Eq. (150) writes in this gauge as 



1 f 2-/V" \ 

d t E PF + d r (r 2 V r £ PF ) = E PF ( — - + d r N r J 



T<)l .(r 2 ^ 3 uVp 



NJ, 



PF 

(r) 



.4 



d r z/ + 3d r a - 3 U {r) AK {r) {r) + NAT 1 



The momentum conservation equations writes 



dt Jr F + N r d r J? F = -J^d r N r -N 



d r PFS [r) [r) + 2 (pfS^) - PF^i)) (- + 5, 
+ (pFS ir) {r) +E PF )d r V+ 3 F r 



n 



In terms of triad components Eq. (|156|) reads, 



PF 

(r) 

JV 
"^4 



iVif ^ - 3V r d r a + iV r (^3d r a 
(£ PF + p) d r v + <9 r p + 3 F( r) A 



d r N r 



In this gauge, the Euler equation (p.57|) reads 



5 S^r) + 3[/ (r) 3o 3 u{ r) = _ 1 [3 a (r) p + S^rjg 

W W -C/PF + P L 

+ ^(.^ ) _. 8tp)1/ ) + _l_(. ff fr)^_ J *)) . 

Under the IGMS coordinates the evolution equations for the entropy per baryon and the 
particle number per baryon keep the same form as Eqs. ( |195[) and ( |197|) or the arternative 
form given by Eqs. fll96D and (li~9gj) , where as Eqs. Q and Q read 



ft (,4 3 n E ) + -i<9 r [r 2 F r n £ 4 
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d t (A 3 n E ) + \d r 



The equation (|226| ) has a conservative form for the quantity n^A? . 

When using the evolution equation (|212j ), the Eqs. ( |226|) and ( |227|) become respectively 



<9 4 n E + r 2 V r n E 



+ 3n E 



iV" 1 

(V r - N r ) 8 r a + -NK r 

r 2 



, 



dtriE + ~^d r 



+ 3h,e 



A 



r 2 



The equation (|22§| ) has a conservative form for the quantity rag. 



B. Radial gauge and polar slicing condition (RGPS) 



In this gauge, B = 1 (/3 = 0) and the polar slicing condition K = K r r is equivalent to 
K e e + K ^ = 0. Since -KT 6 ^ = if ^, this slicing condition and the gauge choice lead to N r = 0. 
The Hamiltonian constraint Eq.( |l34| ) reads 



1 2 

-5 f^ 2 - l) + -d r a = 87rG ^ 2 
V It 



Moreover, by defining 



^M) := 1 



2G m(r,t) 



-1/2 



the Eq.(p30|) reads, 



d r a = A 2 G (inrE - 



or even [cf. Eq.(18) of Ref. || or Eq.(3.29) of Ref. |3(| for a perfect fluid alone or for a perfect 
fluid accompanied by a neutrino flow, respectively], 



d r m = Airr E . 



The momentum constraint (|135| ), for the present gauge choice reads 
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Ri t) = ~^N dtA = 4 7rrG ° Jr = 4 ^ rGo ^ r = 4 ^ rGo A J(r) ' 

This with Eq. ( |231| ) results in an evolution equation for m(r,t) [cf. Eq.(20) of Ref. || or 
Eq.(3.31) of Ref. [|36| for a perfect fluid alone or for a perfect fluid accompanied by a neutrino 
flow, respectively]: 



d t m = -Anr 2 NJ r . 



The evolution equation (|138|) gives, 



d r v 
r 



d r A 
rA 



AttGqA 2 (S (r) ir) -E 



With (EMD and ( Egg) this writes [cf. Eq.(22) of Ref. § and Eq.(3.32) of Ref. [3§ for a source 



term # 2.) of a perfect fluid alone Eq. ( |142|) or that of a perfect fluid accompanied by a 
neutrino flow ( |203j ) , respectively] , 



0,i' : f T ' () A 2 (^ + 4vrr5 (r ( ) r) 



Therefore Eqs. ( |233| ) and (237) are the field equations for the two variables A and N 

respectively. These quantities evolve in time through the matter variables. Therefore, as in 

(r) 

the IGMS coordinate choice, the evolution equation for K A is also redundant in the RGPS 
coordinates. For completeness we write it using Eq. ( |137j ) [cf. Eq.(21) of Ref. H for a perfect 
fluid alone], 



(r) V (r)J j±2 



2d r a 2 



d rr v + (d r u) (d r a — d r v) 



AttG N (-S ir l ) +2S {e) {g) -E 



Concerning the evolution equations for the matter, Eq.( |148D reads, 

a *^PF + ^2 d r ( Ar2 Jp F ) = NK {r) (r) (p F S (r ( } r) + E PF ) - 2J PF d r N - N 2 T f , 



which can be written as 



8 t E PF + ^d r ( Nr 2 J PF ) = NK (r > ( PF S (T > + E PF ) - TV J£ F (d r u + d r a) - 



('■) 



(r) 



T 2-r-t 



(r) 
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We can replace the gradients of the metric potentials by using the Eqs. (E3Q) and ( |236| ) which 
imply, 

2 / o(0 



Then 



d r u + d r a = 4 TrrGo A 2 { S (T > r) + E ) . 



a^pp + (iVr 2 J^ F ) = JVtfWj (p F 5 (r ( } r) + £ PF ) - 4vrrGo JV4 2 J£ F (s^ + ^) - N 2 T l 
Using the momentum constraint (|234| ) , we obtain 
d t E PF + ^d r (Nr 2 J£ F ) = 4vrrG NA 2 [j r ( ppS^ + E PF ) - J r PF (s {r) [r) + E 



N 2j:t 



Now, since 



J — ^PF + , 

= -Epf + -Er , 



* (r) - PF b (r) + (r) 



?(»■) 



7 W 



We find 



d t £ PF + (iVr 2 J? F ) = 4 vrrGo iV"yl 2 [ J R ( ppS^ + £ PF ) - J PF ( rS^ + £r) ] - N 2 T l . 

Using that J PF = ^4Jp F , J R = AJ R and Eqs. ( |142[ ) and ( |144| ) , the energy conservation 
Eq.(p47|) finally reads [cf. Eq.(25) of Ref. § or Eq.(3.55) of Ref. |§ for a perfect fluid alone 
or for a perfect fluid accompanied by a neutrino flow, respectively], 



d t E PF + \d r [r 2 (E PF +p)V r ] = AttvGoNA (E pf + p) [j R r) ((U^) 2 + l) - ( R S {r) {r) + E R 



N 2 P 



The alternative expression of Eq. Q248| ) in conservative form reads 
d t E PF + \d r (r 2 V r E PF ) = -\d r (r 2 V r p)+^rG NA (E PF + p) [j R r) ((U^) 2 + l) - 17« (rS^ + 



The momentum conservation Eq.(|l52j) reads 
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d t J^' + Nd rPF s ir l ) + — PF 



2N 



(r) 
(r) 



pfS 



(9) 
(6) 



PF^ 



S (r l ) +E PF )0,.X+ X •/,''' /\" 



PF 7^(r) 3 



(r) 



This can be written as 



St J r PF = iv 



Using (|156|) the conservative form of this latter reads 



d t J^ )+ U(r^J^)=J { 



PF 

(<') 



2^^ - F r a r a 



X 

1 



(£pf + p) d r u + <9 r p + 3 J"( r )A 



Furthermore, the use of Eqs. (|23l|) , (|3|), (|234|) and (|237| ) lead to 



PF 

(<') 



MrAXJ^ - A 2 V r { 4irrE 



-G NA (E PF +P){^ + ^rS^A - ^ (d rV + 3 F {r) A) . 



m 



N 



Finally, using Eqs. (|244|) and (|246|) and the expressions (|142|) and ([144]), we obtain 



d t 4) + ^d r (rWJf ) = Go Jff ^vrr^ (jff + 2J«) - A 2 V 

-G XA(E FF +p) 



4vrr (£ PF + E R ) 



m 



(r) 



N / , 

-j(drP+ 3 F( r) A 



The Euler equation (157) and Eqs. (234) and ( |237j ) lead to 



d {t) 3 u {r) + 3 c/ (r) 3 d (r) 3 u {r) 



GyA 

r 2 



1 



-EpF + p 



3 d ( - r) p + 3 U {r) d (t)P 



Again, using Eqs. (p44j) and fl246|) and the expressions (142) and (|144|), we obtain 



d t 3 U {r) + V r d r 3 U {r) 



1 '!%, + w>b^ ANa " 



(E FF +p) \A 
+ — - ; ( 3 C/WjpW _ 3^(0 ^ 



r 2 



^+4 7 rr(p +R 5«-3^)j« 



(£ P F + p) 

where we also used that in this gauge = 1/Ndt, dr r \ = 1/Ad r and U^ r > = ^V r . 

This is the Euler equation of the fluid in spherical symmetry which includes the forces of 
the radiation fields acting on the fluid [cf. Eq.(34) of Ref. || or Eq.(3.56) of Ref. Q| for a 
perfect fluid alone or for a perfect fluid accompanied by a neutrino flow, respectively]. 
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A relation that turns to be useful in this gauge is obtained by combining Eqs. (234) and 



m 



AN 



d t (A 2 ) + d r v + d r a = A-ktGqA 2 



S {r) [r) + E- 2 JM 



For example, in the static case and for perfect fluids, this provides a simple relation between 
the gradients of the metric potentials and the pressure and energy-density of matter. Moreover, 
outside the star surface the only contributions to the total matter variabels are those of radiated 
particles. Under the free streaming approximation [cf. Eq.(187)], this implies that the r.h.s of 
fl257|) vanishes. This situation was investigated analitically in the past using a different gauge 



|69|j70| and corresponds to the external solution. 

Finally, the integrability condition d 2 , t m = d 2 r m imposed in Eqs. (|233| ) and (|235| ) result in 
the relationship 



1 



d t E + -2 d r [r 2 NJ r j = . 

This equation is in fact compatible with the evolution equation for the total energy density of 
matter. 

Indeed, subtracting Eq. (|243| ) or more specifically Eq. fl248|) from (|258 ) one obtains an 
evolution equation for the energy-density of radiation Er [cf . Eq. (|191| ) ] . 

The evolution equation ( |113| ) and the alternative form ( |194J ) in RGPS coordinates read 
respectively 

d t (An E ) + \d r [rVAn E ] = , 



d t (An B ) + \d r (r 2 n E N 3 U^) = 



Note that Eq. Q25S| ) has a conservative form for the quantity An E . 
Using Eqs. pi), (|232|), (J24§), (|U3) and p6|) in Eq. (ggg) it turns 



r 2 n E V r 



GqA TIE 



m 



V r + Aixrp + AirrNJ] 



R 
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which provides an equation in conservative form for n E . 



d t n E + -^d r (r 2 3 U^n E ) + n E G Q AN 



In the same way, using Eqs. ( |234| ) and ( |237|) in the alternative Eq. ( 260 ) we obtain 

= (262) 



The evolution Eqs. ( |195D and ( |197| ) do not change in form under the RGPS coordinates. 
However, the alternative form given by Eqs. ( |196| ) and ( |198| ) in RGPS coordinates write 
respectively 

d t a+ d r a = ——U n n R + V p ) . (263) 

A nl B V / 

AT 3 fjW N 

dtXR + d r x R = —K R . (264) 

A nl 



C. Isotropic gauge and polar slicing condition (IGPS) 

The isotropic choice A = B (i.e., a = f3), and the polar slicing condition K e e + K^^ = 0, 
implies due to the spherical symmetry that K e g = 0. This latter leads to an evolution equation 
for a: 

N r 

d t a = -N r d r a . (265) 

r 



This and the expression for K r r [cf. Eq, (|12S| )l provide an equation for the shift: 



/ - / / - v) 



d r{ t2 -)= -^^t) ■ ( 2 66) 



The Hamiltonian constraint (|134j ) leads to 

5 r 2 r a + — — = -2vrGo^ 2 ^ - (9 r a) 2 , (267) 
r 

where, 

q = q/2 . (268) 
The momentum constraint (|135| ) writes, 
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('') 



1 



K K L ( - + d r a I = 4 7rG A J (r) = 4vrG J r • 



The equation (|13§| ) together with fl267|) provide an equation for the lapse, 



1 



i),v [ - + 2f; ( .ri ) = 4?r(?oA 2 5 (r ( J. ) - 20 r a | - + r>,.o 



In the case of a perfect fluid alone, this corresponds to Eq. (9) of Shapiro & Teukolsky 4C], 
and to Eq.(6) of Schinder et al. where the authors use a "Lagrangian" gauge which can be 
easily transformed to the IGPS coordinates. 

Finally, the evolution equation for K ) * given by Eq. (|137|) together with Eq. fl267|) yield 



9tK% + 



(r) 

, Mf N 



-d TT v + [d r v) (d r a — d r v) + (d r a) — h 2d r a 



= 4ttG N (-S (r l ) +2S i % ) +E) . 
The equation of conservation of energy ( |150|) reads in this gauge as follows, 
Q tEpF + ^ dr (r 2 ^ E PF ) = E PF (NK^ r) + ^ + d r N^ - ±d r (r 2 ^ 3 [/Wp) 



PF 

(r) 



<V + 3d r a - 3 U {r) AK (r) (r) 



A 



The equation of conservation of momentum Eq. (|156| ) reads 



+ NK ir) (r) p - N 2 T l . 



PF 

(r) 

iV 
'A 



2NK {r ^ r) - W r d r a + +N r f :W,.n - - ) - ( ';,.Y 



(-EpF + p) <V + <9 r p + 3 .F( r )A 



The Euler equation (157) in this gauge take the same form of equation (|225j), and the 



Eqs. (|196|), (|198|) and (194) [see Eqs. (|226|) and ([227])] keep also the same form as in the 



IGMS coordinates. Altenatively, one can also use the simpler form of Eqs. ( |195| ) and (197). 
Furthermore, when using the evolution equation ( |265| ) in Eqs. (199) and ( |199| ) in IGPS 
coordinates we obtain respectively 



1 



dtfiE + -*d r r 2 V r nE + 3ng 



jV r ' 

fyr _ N rj drQL 

r 



, 



d t n E + -^d r 



+ 3ue 



A 



. 
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D. Boundary conditions and initial data 



A typical feature of spherically symmetric spacetimes is that the gravitational field variables 
can evolve in time through the matter fields. So the initial conditions for the matter variables 
and the boundary conditions fix automatically the initial values for the gravitational field 
by solving the constraint equations and the slicing condition equation. For spacetimes with 
less symmetries one is always forced to solve the dynamic Einstein equations to evolve the 
gravitational field. Let us thus discuss first the boundary conditions. 

We call exterior solution the solution of field equations outside the perfect-fluid domain 
(usually a compact support). In the present case, it does not correspond to the Schwarzschild 
vacuum solution since in general, the radiated matter will extend to spatial infinity. Thus, the 
exterior solution has to be found also numerically. The exterior sources of the field equations 



will be provided by the energy-momentum tensor of particles [cf. Eq. (176)]. The matter 
variables of particles will evolve in time through the distribution function. Moreover, outside 
the star, the radiated particles can interact only with themselves, however, this interaction is 
rather week in comparison with the interaction inside the star. In a first approximation one 
can thus neglect such interactions and consider that the particles will follow geodesies; the 
distribution function will thus remain constant along them. 

Regarding the boundary conditions, these are rather regularity and asymptotic conditions. 
For instance, the regularity and the asymptotic flatness condition for the shift are respectively 
(see Ref. [63] for a more detailed analysis about regularity and boundary conditions) 

N r (t,0) = , (276) 
iV r (t,r)™^ . (277) 



Similar conditions apply for K r r . These boundary conditions are enforced from Eqs. fl211| ) 
and ( |266| ) : 

N r (t, r) = A TS r J™ (t, r>)dr> (278) 
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where Ats = 1)3/2 for the IGPS and IGMS coordinates respectively. 

The condition for the lapse at the star's center is such that the asymptotic flatness condition 
N — > 1 is verified. Therefore 

N(t,0) = N c (t) , (279) 

with N c (t) such that 

iV^r)™ ^ 1 . (280) 

Since apriori this is difficult to enforce, a better strategy consists in rescaling the lapse as 
N = N/N c so that N c (t) = 1, then the values of the true lapse can be recovered by using 
N c (t) = 1/N(t, r) r _ >00 where the asymptotic value Noo(t) := N(t, r),.-^ is found numerically 
at every time step. This rescaling allows to integrate the equations spatially in only one 
cycle. The rescaling will not affect the relevant equations of motion for N or A since only the 
derivatives of v appear there [i.e., the Eqs. pTBD , ( gig ), (ggg), fl236| ), (|267| ), ( p70|) for iV and 



A in the different gauges are invariant to such a rescaling]. However, this is not true for the 
shift equation and for the Boltzmann equation where N appears explicitly. However, this does 
not pose a problem since a simultaneous rescaling N r = N r /N c leaves all equations invariant 
as well as the boundary conditions for N r . 

In the case of the RGPS and IGPS coordinates one can find an integral expression for the 
lapse satisfying the boundary conditions. For instance, from Eq. fl237p 



u(t, r) = G £ A 2 + 4tt r' dr' + u(t, 0) . (281) 

The asymptotic flatness condition Eq. (|280| ) leads to 

v(t,r) r -+oc -» . (282) 
Therefore from Eq. ( |281| ) and the asymptotic condition one obtains 

u(t, 0) = -Go J q °° A 2 + An r' S {r) {r ^ dr' . (283) 
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This corresponds precisely to the renormalized value — u^. So finally, 

u(t, r) = -Go jH A 2 + 4tt r' S {r) {r ^ dr' . (284) 

The value for the lapse at the star surface will be provided by 

u(t, R(t)) = -G jH A 2 + 4vr r' dr' , (285) 

where R(t) corresponds to the RGPS-r coordinate at the star surface at time t. We emphasize 
that, outside the star S^ r K = rS 1 ^, that is, the only contribution to S vL is from the radiated 
particles. In fact, outside the star we can write, 

v&rU = -Go£ m A 2 {^+^r' R S^ dr' 

= []nA]% m - 4vr A 2 r' + E R ) dr' 

= lnA^r)- 1 !^) - 4tt jT* ^ A 2 r' ( R 5 (r ( } r) + E R ) dr' 



In 



1 



2G m(t,r^ 1/2 



r>R(t) 

4tt A 2 r' ( K S^L + E R ) dr' , (286) 

Jr>R(t) V ( > ' 



r J r>R(i) Jr>R(t) 
where we used Eq.( p32| ) in order to replace the term with m/r' 2 and also the asymptotic 
flatness condition on A [cf. Eq.( |289| ) below]. It is to be stressed that in the absence of matter 
outside the star surface, the first term of Eq.( ^86|) corresponds to the expression for v of the 
Schwarzshild metric (with m(t, t) | r >ii(t) = being the total mass of the star). However, 
in the present case there are contributions (due to the energy-density E R and pressure rS^^ 
of particles which fills some part of the space outside the star) which are the responsible for 
the actual metric to deviate from the Schwarzschild one. These contributions arise from the 



second term of Eq. (286). In some cases (e.g. the free streaming approximation) numerical 
analysis show that these deviations are so small so that they can be neglected (cf . [|jfj ) . 

Deviations of this kind can be appreciated more easily in presence of non trivial scalar fields, 
for instance in the phenomenon of spontaneous scalarization [71|. Moreover, for r > R(t), the 
mass function m(t, r) is larger than m(t, R(t)) due to the contribution of E R to the total energy 
density [cf. Eq. (|233|) ]. Indeed the mass difference is given by 
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rr>R(t) 

5m = A-K / E R r' 2 dr' . (287) 

JR(t) 

Another way to appreciate such deviations is by noting that Eq. ( |257| ) together with 
the asymptotic conditions imply, in the case of vacuum, the relationship AN = 1 which is 
characteristic of the Schwarzschild solution in RGPS coordinates. However, when matter is 
present outside the star, then AN ^ 1 (e.g. see ffljl), except of course at spatial infinity 

The boundary conditions for A are similar to those for N. Therefore 

A{t,0) = A c (t) . (288) 

with A c (t) such that 

A(i,r) r ^oo -> 1 . (289) 



In the case of the RGPS coordinates the reparametrization Eq. (|231|) imposes the regularity 
condition 

m(t,0) = . (290) 

Since near the origin m ~ r 3 , the reparametrization inforces that A(t, 0) = 1. Then dtA(t, 0) = 
0. The three metric is thus locally flat at the origin. Moreover, provided that the energy-density 
of sources falls off at least as fast as 1/r outside the star, the same mass parametrization will 
drive A to the required asymptotic value. Note, that this behavior of the metric potential A is 
compatible with the regularity and asymptotic conditions for K r r that we mentioned above. 
Moreover, these conditions also imply that J r vanish at the origin as well as asymptotically 
[cf. Eq.([23D ]. 

In the case of the IGMS and IGPS coordinates one has an elliptic equation for a which is 
not invariant to a rescaling on A. So unlike the RGPS coordinates where A c = 1, the central 
value A c , should be determined from a shooting method or otherwise in order to satisfy the 
asymptotic flatness condition. In fact, near the origin N r ~ r, therefore from Eq. ( |212| ), it 
turns that at the origin 
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d t a(t, 0) = const. , (291) 
and thus from the definition Eq. ( |132| ) and the regularity condition ( |292| ) (see below) one con- 



cludes that (depending on the sign of the constant) A(t, 0) can grow exponentially flq jlC||3q{l2 | . 
The three metric is thus conformally flat at the origin in the isotropic gauge. 
In addition to the previous regularity conditions we have also 

a r Q| r=0 = = Q r |r=o • (292) 

where Q represents collectively the metric potentials and the scalar matter field variables (e.g, 
N, A,m,p, p, E, etc.) and Q r tensor field components (e.g J r ,K r r , etc.). 

A convenient way to impose the asymptotic conditions accurately is by compactifying the 
outer domain with the help of a transformation u = 1/r from the star's surface r = R to 
spatial infinity. In this way the infinite domain r £ [R, oo) is mapped to the compact domain 
u € [1/i?, 0), so the integration can be performed from 1/R to "zero" with a high degree of 
accuracy (cf. [[n]]). Obviously, the field and matter variables for particles are to be matched 
continously at R. 

Regarding the distribution function, the regularity condition at the center on the particle's 
radial energy flux is H R = [see Eq. ([1781 )1 which means that the average of the particles' 



radial velocity as measured by the Eulerian observer at the origin is zero (local isotropy at 
r = 0). Same considerations apply for the radial particle-flux- number j( r ) [cf. Eq. ( |190| )1. A 
sufficient condition for H R = = at r = is 

F R {t, 0, e, //) = G R {t, e, p) with G R {t, e, y) = G R (t, e, -//) , (293) 

that is, the distribution function being a pair function of \x at r = 0, enforces that the integrals 
given by Eqs. ( |17§| ) and ( |190| ) vanish identically at the origin. This condition ensures in 
addition that with a suitable form of G R (t,e,n) in the energy domain the energy-density, 
pressure and density number of particles in the Eulerian frame are finite at the star's center. 
Usually the assumption of an ideal quantum gas is adopted as initial condition for the particles 
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so that the distribution function is isotropic and given by a Fermi-Dirac or Bose-Einstein 
distribution (whether the particles are Fermions or Bosons) |p0|,^ , 45 ] . Then the particles 



energy-density and pressure will be parametrized initially only by the thermodynamic variables 
like the temperature. In that case Hr = all over the initial spacelike hypersurface. 
The regularity condition for F R at r = is like other scalar quantities, 

d r F R \ r=0 = . (294) 

Another boundary condition is that the inward flux of particles at the star surface is zero. 
This is imposed by |56| , 

F R \ r=R = for - 1 < n < . (295) 

Returning to the initial conditions, the form of these, will characterize first, the type of configu- 
ration (e.g., supernova, neutron star, supermassive star, star cluster) and second the dynamics 
and ultimate fate of the precursor (e.g. protoneutron star, neutron star, black hole). The goal 
of future numerical investigation will be to explore a large set of initial conditions and their 
consequences (cf. fj%). 



VIII. CONCLUSIONS 

Several astrophysical phenomena involve the dynamics of relativistic objects. Some of the 
most interesting ones end up in the formation of black holes or neutron stars, like the collapse 
of cores and supernova explosions. While most of the astrophysical objects are rotating, the 
role of rotation in relativity can be neglected as regards the structure of the object when 
the rotation frequency is low. Aside from fast pulsars, most of the astrophysical objects are 
endowed with a low rotation frequency. Therefore, the spherical symmetry is an assumption 
that can be very useful in a wide range of applications. On the other hand, the mere existence 
of fast pulsars reveals that rotation has to be taken into acount in a more realistic situation. 
Moreover, it seems that the deviations from spherical symmetry in supernova explosions is 
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central in the phenomenon jlT]], and that rotation can influence the cooling mechanims in 



early phases of neutron stars [73|. In general relativity this is a difficult task and only a few 



attempts have been succeeded within an evolutionary code (see [74 and references therein). 

One can separate the problem of the dynamics of relativistic bodies in two sets. The 
first one involves the formulation, the geometry and the numerics. A convenient coordinate 
choice and a powerful numerical analysis can allow long term evolutions leading to a better 
understanding of several physical phenomena. Thus this is a crucial point which has been 
recognized by almost all the numerical relativists. Investigations of the effect of several gauges 
and time slicing conditions is always an important issue. One of the aims of this paper was 
therefore to derive the fundamental equations under different gauges and write them in several 
forms suitable for different numerical schemes. 

The second set, involves the physical approximations used in the model. In the case of 
gravitational collapse, we have discussed that neutrinos cannot only play an important role 
in the dynamics but also that the signal carried by them can be fundamental for a better 
understanding of the underlying physics and as an invaluable imprint of the ultimate fate of 



the collapsed object. In particular, if neutrinos turn to be massive particles [34], mechanisms 



like the early black hole formation could be tested [22|. 

Furthermore, the equation of state at high densities can also lead to different time scale 
processes during the collapse and the accretion phase. Therefore, it turns necessary to pursue 
the analysis with the incorporation of the most recent advances in particle and nuclear physics. 

While the formalism presented here included only hydrodynamics and transport theory, 
the equations are quiet general as to include other types of matter like scalar fields, which are 
very useful in the analysis of critical phenomena. 

Our aim for the future investigations is to perform an extensive numerical analysis of several 
issues discussed here and more generally to analyse the dynamics of spherically symmetric 
spacetimes with several kind of matter sources. 
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